function [eqmin,FeMgKd] = fractionateOlivine_V2(liquidComp,FeMgKd)

%FeMgKd = 0.29;  % Kd = (Fe/Mg)min / (Fe/Mg)liq
 
 % Columns are Ti, Cr, K, Na
% Rows are cpx, opx, olv, plagioclase/spinel/garnet
 
% DMP = ...
%     [0.42	4.3       0.001	    0.15 
%     0.19	5         0.001	    0.04
%     0.03	1.45     0.001	   0.001
%     0.001   0.001	0.271	  1]
 
 D_TiO2 = 0.03;
 D_Cr2O3 = 1.45;
 D_K2O = 0.001;
 D_Na2O = 0.001;
 
MW = [60.085 79.899	101.961	151.9902	71.846	70.9374	40.311	56.079	61.979	94.203	283.89]; % see matlab guide in excel spreadsheet
%     1-SiO2    2-TiO2    3-Al2O3   4-Cr2O3       5-FeO     6-MnO     7-MgO     8-CaO    9-Na2O    10-K2O    P2O5

% use molar weights and Kds to calculate ratios of elements in mineral

liquidComp_moles = liquidComp./MW; 
liquidComp_moles = liquidComp_moles/nansum(liquidComp_moles).*100;


%   W = 3000;
%   R = 8.31446;
%   SiO2_Toplis = liquidComp_moles(1);
%   XFo_Toplis = .88;
%   T = 1200; 
%   P = 10; 
%   
%   
%   T_Toplis = T;
%   P_Toplis = P.*10^3;
%   
%   
%   KD_Toplis = exp((-6766./R./T_Toplis - 7.34./R) + log(0.036.*SiO2_Toplis-0.22)+(W./R./T_Toplis.*(1-2.*XFo_Toplis))+...
%       (0.035./R./T_Toplis.*(P_Toplis-1))) %% Toplis KD prediction
%   
  
  
  XSi = liquidComp_moles(1); 
  XNaK = liquidComp_moles(9)+liquidComp_moles(10); 
%   XSi./100
%   XNaK./100
  Longhi_KD = -.0804 + 0.542.*XSi./100 - 0.594.*XNaK./100;
  Longhi_KD2 = 0.165 + 0.299.*XSi./100 ;
  
  FeMgKd = Longhi_KD2; 
  %FeMgKd = Longhi_KD; 
  
    FeMgliq = (liquidComp(5)./MW(5))./ (liquidComp(7)./MW(7)); %% could instead first convert liquid to moles, then calculate the mineral comp. 
    FeMgmin = FeMgKd*FeMgliq;
    Xol = FeMgmin/(1+FeMgmin);
    

    M_eqmin(1) = 33.33;
    M_eqmin(2) = 0;
    M_eqmin(3) = 0;
    M_eqmin(4) = 0;
    M_eqmin(5) = 66.67*Xol;
    M_eqmin(6) = 0;
    M_eqmin(7) = 66.67*(1-Xol);
    M_eqmin(8) = 0;
    M_eqmin(9) = 0;
    M_eqmin(10) = 0;
    M_eqmin(11) = 0;
    


    eqmin = M_eqmin.*MW;                % this is the bulk comosition of equilibrium olivine
    eqmin = eqmin./(0.01*nansum(eqmin));   % renormalized and in wt%
    
    
    eqmin(2) = D_TiO2.*liquidComp(2);
    eqmin(4) = D_Cr2O3.*liquidComp(4);
    eqmin(9) = D_Na2O.*liquidComp(9);
    eqmin(10) = D_K2O.*liquidComp(10);
    eqmin = eqmin./(0.01*nansum(eqmin));   % renormalized and in wt%

end

